/* linalg/gsl_linalg.h
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_LINALG_FLOAT_H__
#define __GSL_LINALG_FLOAT_H__

#include <gsl/gsl_mode.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_vector_float.h>
#include <gsl/gsl_matrix_float.h>
#include <gsl/gsl_linalg.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
#define __BEGIN_DECLS extern "C" {
#define __END_DECLS }
#else
#define __BEGIN_DECLS           /* empty */
#define __END_DECLS             /* empty */
#endif

__BEGIN_DECLS

/* typedef enum */
/*   { */
/*     GSL_LINALG_MOD_NONE = 0, */
/*     GSL_LINALG_MOD_TRANSPOSE = 1, */
/*     GSL_LINALG_MOD_CONJUGATE = 2 */
/*   } */
/* gsl_linalg_matrix_mod_t; */


/* Note: You can now use the gsl_blas_dgemm function instead of matmult */

/* Simple implementation of matrix multiply.
 * Calculates C = A.B
 *
 * exceptions: GSL_EBADLEN
 */
int gsl_linalg_matmult_float (const gsl_matrix_float * A,
                              const gsl_matrix_float * B,
                              gsl_matrix_float * C);


/* Simple implementation of matrix multiply.
 * Allows transposition of either matrix, so it
 * can compute A.B or Trans(A).B or A.Trans(B) or Trans(A).Trans(B)
 *
 * exceptions: GSL_EBADLEN
 */
int gsl_linalg_matmult_mod_float (const gsl_matrix_float * A,
                                  gsl_linalg_matrix_mod_t modA,
                                  const gsl_matrix_float * B,
                                  gsl_linalg_matrix_mod_t modB,
                                  gsl_matrix_float * C);

/* Calculate the matrix exponential by the scaling and
 * squaring method described in Moler + Van Loan,
 * SIAM Rev 20, 801 (1978). The mode argument allows
 * choosing an optimal strategy, from the table
 * given in the paper, for a given precision.
 *
 * exceptions: GSL_ENOTSQR, GSL_EBADLEN
 */
int gsl_linalg_exponential_ss_float (
                                     const gsl_matrix_float * A,
                                     gsl_matrix_float * eA,
                                     gsl_mode_t mode
                                     );


/* Householder Transformations */

double gsl_linalg_householder_transform_float (gsl_vector_float * v);

int gsl_linalg_householder_hm_float (double tau, 
                                     const gsl_vector_float * v, 
                                     gsl_matrix_float * A);

int gsl_linalg_householder_mh_float (double tau, 
                                     const gsl_vector_float * v, 
                                     gsl_matrix_float * A);

int gsl_linalg_householder_hv_float (double tau, 
                                     const gsl_vector_float * v, 
                                     gsl_vector_float * w);

int gsl_linalg_householder_hm1_float (double tau, 
                                      gsl_matrix_float * A);

/* Singular Value Decomposition

* exceptions: 
*/

int
gsl_linalg_SV_decomp_float (gsl_matrix_float * A,
                            gsl_matrix_float * V,
                            gsl_vector_float * S,
                            gsl_vector_float * work);

int
gsl_linalg_SV_decomp_mod_float (gsl_matrix_float * A,
                                gsl_matrix_float * X,
                                gsl_matrix_float * V,
                                gsl_vector_float * S,
                                gsl_vector_float * work);

int gsl_linalg_SV_decomp_jacobi_float (gsl_matrix_float * A,
                                       gsl_matrix_float * Q,
                                       gsl_vector_float * S);

int
gsl_linalg_SV_solve_float (const gsl_matrix_float * U,
                           const gsl_matrix_float * Q,
                           const gsl_vector_float * S,
                           const gsl_vector_float * b,
                           gsl_vector_float * x);


/* LU Decomposition, Gaussian elimination with partial pivoting
 */

int gsl_linalg_LU_decomp_float (gsl_matrix_float * A,
                                gsl_permutation * p,
                                int *signum);

int gsl_linalg_LU_solve_float (const gsl_matrix_float * LU,
                               const gsl_permutation * p,
                               const gsl_vector_float * b,
                               gsl_vector_float * x);

int gsl_linalg_LU_svx_float (const gsl_matrix_float * LU,
                             const gsl_permutation * p,
                             gsl_vector_float * x);

int gsl_linalg_LU_refine_float (const gsl_matrix_float * A,
                                const gsl_matrix_float * LU,
                                const gsl_permutation * p,
                                const gsl_vector_float * b,
                                gsl_vector_float * x,
                                gsl_vector_float * residual);

int gsl_linalg_LU_invert_float (const gsl_matrix_float * LU,
                                const gsl_permutation * p,
                                gsl_matrix_float * inverse);

double gsl_linalg_LU_det_float (gsl_matrix_float * LU, int signum);
double gsl_linalg_LU_lndet_float (gsl_matrix_float * LU);
int gsl_linalg_LU_sgndet_float (gsl_matrix_float * lu, int signum);

/* QR decomposition */

int gsl_linalg_QR_decomp_float (gsl_matrix_float * A,
                                gsl_vector_float * tau);

int gsl_linalg_QR_solve_float (const gsl_matrix_float * QR,
                               const gsl_vector_float * tau,
                               const gsl_vector_float * b,
                               gsl_vector_float * x);

int gsl_linalg_QR_svx_float (const gsl_matrix_float * QR,
                             const gsl_vector_float * tau,
                             gsl_vector_float * x);

int gsl_linalg_QR_lssolve_float (const gsl_matrix_float * QR, 
                                 const gsl_vector_float * tau, 
                                 const gsl_vector_float * b, 
                                 gsl_vector_float * x, 
                                 gsl_vector_float * residual);


int gsl_linalg_QR_QRsolve_float (gsl_matrix_float * Q,
                                 gsl_matrix_float * R,
                                 const gsl_vector_float * b,
                                 gsl_vector_float * x);

int gsl_linalg_QR_Rsolve_float (const gsl_matrix_float * QR,
                                const gsl_vector_float * b,
                                gsl_vector_float * x);

int gsl_linalg_QR_Rsvx_float (const gsl_matrix_float * QR,
                              gsl_vector_float * x);

int gsl_linalg_QR_update_float (gsl_matrix_float * Q,
                                gsl_matrix_float * R,
                                gsl_vector_float * w,
                                const gsl_vector_float * v);

int gsl_linalg_QR_QTvec_float (const gsl_matrix_float * QR,
                               const gsl_vector_float * tau,
                               gsl_vector_float * v);

int gsl_linalg_QR_Qvec_float (const gsl_matrix_float * QR,
                              const gsl_vector_float * tau,
                              gsl_vector_float * v);

int gsl_linalg_QR_unpack_float (const gsl_matrix_float * QR,
                                const gsl_vector_float * tau,
                                gsl_matrix_float * Q,
                                gsl_matrix_float * R);

int gsl_linalg_R_solve_float (const gsl_matrix_float * R,
                              const gsl_vector_float * b,
                              gsl_vector_float * x);

int gsl_linalg_R_svx_float (const gsl_matrix_float * R,
                            gsl_vector_float * x);


/* Q R P^T decomposition */

int gsl_linalg_QRPT_decomp_float (gsl_matrix_float * A,
                                  gsl_vector_float * tau,
                                  gsl_permutation * p,
                                  int *signum,
                                  gsl_vector_float * norm);

int gsl_linalg_QRPT_decomp2_float (const gsl_matrix_float * A, 
                                   gsl_matrix_float * q, gsl_matrix_float * r, 
                                   gsl_vector_float * tau, 
                                   gsl_permutation * p, 
                                   int *signum,
                                   gsl_vector_float * norm);

int gsl_linalg_QRPT_solve_float (const gsl_matrix_float * QR,
                                 const gsl_vector_float * tau,
                                 const gsl_permutation * p,
                                 const gsl_vector_float * b,
                                 gsl_vector_float * x);


int gsl_linalg_QRPT_svx_float (const gsl_matrix_float * QR,
                               const gsl_vector_float * tau,
                               const gsl_permutation * p,
                               gsl_vector_float * x);

int gsl_linalg_QRPT_QRsolve_float (const gsl_matrix_float * Q,
                                   const gsl_matrix_float * R,
                                   const gsl_permutation * p,
                                   const gsl_vector_float * b,
                                   gsl_vector_float * x);

int gsl_linalg_QRPT_Rsolve_float (const gsl_matrix_float * QR,
                                  const gsl_permutation * p,
                                  const gsl_vector_float * b,
                                  gsl_vector_float * x);

int gsl_linalg_QRPT_Rsvx_float (const gsl_matrix_float * QR,
                                const gsl_permutation * p,
                                gsl_vector_float * x);

int gsl_linalg_QRPT_update_float (gsl_matrix_float * Q,
                                  gsl_matrix_float * R,
                                  const gsl_permutation * p,
                                  gsl_vector_float * u,
                                  const gsl_vector_float * v);

/* LQ decomposition */

int gsl_linalg_LQ_decomp_float (gsl_matrix_float * A, gsl_vector_float * tau);

int gsl_linalg_LQ_solve_T_float (const gsl_matrix_float * LQ, const gsl_vector_float * tau, 
                                 const gsl_vector_float * b, gsl_vector_float * x);

int gsl_linalg_LQ_svx_T_float (const gsl_matrix_float * LQ, const gsl_vector_float * tau, 
                               gsl_vector_float * x);

int gsl_linalg_LQ_lssolve_T_float (const gsl_matrix_float * LQ, const gsl_vector_float * tau, 
                                   const gsl_vector_float * b, gsl_vector_float * x, 
                                   gsl_vector_float * residual);

int gsl_linalg_LQ_Lsolve_T_float (const gsl_matrix_float * LQ, const gsl_vector_float * b, 
                                  gsl_vector_float * x);

int gsl_linalg_LQ_Lsvx_T_float (const gsl_matrix_float * LQ, gsl_vector_float * x);

int gsl_linalg_L_solve_T_float (const gsl_matrix_float * L, const gsl_vector_float * b, 
                                gsl_vector_float * x);

int gsl_linalg_LQ_vecQ_float (const gsl_matrix_float * LQ, const gsl_vector_float * tau, 
                              gsl_vector_float * v);

int gsl_linalg_LQ_vecQT_float (const gsl_matrix_float * LQ, const gsl_vector_float * tau, 
                               gsl_vector_float * v);

int gsl_linalg_LQ_unpack_float (const gsl_matrix_float * LQ, const gsl_vector_float * tau, 
                                gsl_matrix_float * Q, gsl_matrix_float * L);

int gsl_linalg_LQ_update_float (gsl_matrix_float * Q, gsl_matrix_float * R,
                                const gsl_vector_float * v, gsl_vector_float * w);
int gsl_linalg_LQ_LQsolve_float (gsl_matrix_float * Q, gsl_matrix_float * L, 
                                 const gsl_vector_float * b, gsl_vector_float * x);

/* P^T L Q decomposition */

int gsl_linalg_PTLQ_decomp_float (gsl_matrix_float * A, gsl_vector_float * tau, 
                                  gsl_permutation * p, int *signum, 
                                  gsl_vector_float * norm);

int gsl_linalg_PTLQ_decomp2_float (const gsl_matrix_float * A, gsl_matrix_float * q, 
                                   gsl_matrix_float * r, gsl_vector_float * tau, 
                                   gsl_permutation * p, int *signum, 
                                   gsl_vector_float * norm);

int gsl_linalg_PTLQ_solve_T_float (const gsl_matrix_float * QR,
                                   const gsl_vector_float * tau,
                                   const gsl_permutation * p,
                                   const gsl_vector_float * b,
                                   gsl_vector_float * x);

int gsl_linalg_PTLQ_svx_T_float (const gsl_matrix_float * LQ,
                                 const gsl_vector_float * tau,
                                 const gsl_permutation * p,
                                 gsl_vector_float * x);

int gsl_linalg_PTLQ_LQsolve_T_float (const gsl_matrix_float * Q, const gsl_matrix_float * L,
                                     const gsl_permutation * p,
                                     const gsl_vector_float * b,
                                     gsl_vector_float * x);

int gsl_linalg_PTLQ_Lsolve_T_float (const gsl_matrix_float * LQ,
                                    const gsl_permutation * p,
                                    const gsl_vector_float * b,
                                    gsl_vector_float * x);

int gsl_linalg_PTLQ_Lsvx_T_float (const gsl_matrix_float * LQ,
                                  const gsl_permutation * p,
                                  gsl_vector_float * x);

int gsl_linalg_PTLQ_update_float (gsl_matrix_float * Q, gsl_matrix_float * L,
                                  const gsl_permutation * p,
                                  const gsl_vector_float * v, gsl_vector_float * w);

/* Cholesky Decomposition */

int gsl_linalg_cholesky_decomp_float (gsl_matrix_float * A);

int gsl_linalg_cholesky_solve_float (const gsl_matrix_float * cholesky,
                                     const gsl_vector_float * b,
                                     gsl_vector_float * x);

int gsl_linalg_cholesky_svx_float (const gsl_matrix_float * cholesky,
                                   gsl_vector_float * x);


/* Cholesky decomposition with unit-diagonal triangular parts.
 *   A = L D L^T, where diag(L) =_float (1,1,...,1).
 *   Upon exit, A contains L and L^T as for Cholesky, and
 *   the diagonal of A is_float (1,1,...,1). The vector Dis set
 *   to the diagonal elements of the diagonal matrix D.
 */
int gsl_linalg_cholesky_decomp_unit_float(gsl_matrix_float * A, gsl_vector_float * D);


/* Symmetric to symmetric tridiagonal decomposition */

int gsl_linalg_symmtd_decomp_float (gsl_matrix_float * A, 
                                    gsl_vector_float * tau);

int gsl_linalg_symmtd_unpack_float (const gsl_matrix_float * A, 
                                    const gsl_vector_float * tau,
                                    gsl_matrix_float * Q, 
                                    gsl_vector_float * diag, 
                                    gsl_vector_float * subdiag);

int gsl_linalg_symmtd_unpack_T_float (const gsl_matrix_float * A,
                                      gsl_vector_float * diag, 
                                      gsl_vector_float * subdiag);

/* Linear Solve Using Householder Transformations

* exceptions: 
*/

int gsl_linalg_HH_solve_float (gsl_matrix_float * A, const gsl_vector_float * b, gsl_vector_float * x);
int gsl_linalg_HH_svx_float (gsl_matrix_float * A, gsl_vector_float * x);

/* Linear solve for a symmetric tridiagonal system.

* The input vectors represent the NxN matrix as follows:
*
*     diag[0]  offdiag[0]             0    ...
*  offdiag[0]     diag[1]    offdiag[1]    ...
*           0  offdiag[1]       diag[2]    ...
*           0           0    offdiag[2]    ...
*         ...         ...           ...    ...
*/
int gsl_linalg_solve_symm_tridiag_float (const gsl_vector_float * diag,
                                         const gsl_vector_float * offdiag,
                                         const gsl_vector_float * b,
                                         gsl_vector_float * x);

/* Linear solve for a nonsymmetric tridiagonal system.

* The input vectors represent the NxN matrix as follows:
*
*       diag[0]  abovediag[0]              0    ...
*  belowdiag[0]       diag[1]   abovediag[1]    ...
*             0  belowdiag[1]        diag[2]    ...
*             0             0   belowdiag[2]    ...
*           ...           ...            ...    ...
*/
int gsl_linalg_solve_tridiag_float (const gsl_vector_float * diag,
                                    const gsl_vector_float * abovediag,
                                    const gsl_vector_float * belowdiag,
                                    const gsl_vector_float * b,
                                    gsl_vector_float * x);


/* Linear solve for a symmetric cyclic tridiagonal system.

* The input vectors represent the NxN matrix as follows:
*
*      diag[0]  offdiag[0]             0   .....  offdiag[N-1]
*   offdiag[0]     diag[1]    offdiag[1]   .....
*            0  offdiag[1]       diag[2]   .....
*            0           0    offdiag[2]   .....
*          ...         ...
* offdiag[N-1]         ...
*/
int gsl_linalg_solve_symm_cyc_tridiag_float (const gsl_vector_float * diag,
                                             const gsl_vector_float * offdiag,
                                             const gsl_vector_float * b,
                                             gsl_vector_float * x);

/* Linear solve for a nonsymmetric cyclic tridiagonal system.

* The input vectors represent the NxN matrix as follows:
*
*        diag[0]  abovediag[0]             0   .....  belowdiag[N-1]
*   belowdiag[0]       diag[1]  abovediag[1]   .....
*              0  belowdiag[1]       diag[2]
*              0             0  belowdiag[2]   .....
*            ...           ...
* abovediag[N-1]           ...
*/
int gsl_linalg_solve_cyc_tridiag_float (const gsl_vector_float * diag,
                                        const gsl_vector_float * abovediag,
                                        const gsl_vector_float * belowdiag,
                                        const gsl_vector_float * b,
                                        gsl_vector_float * x);


/* Bidiagonal decomposition */

int gsl_linalg_bidiag_decomp_float (gsl_matrix_float * A, 
                                    gsl_vector_float * tau_U, 
                                    gsl_vector_float * tau_V);

int gsl_linalg_bidiag_unpack_float (const gsl_matrix_float * A, 
                                    const gsl_vector_float * tau_U, 
                                    gsl_matrix_float * U, 
                                    const gsl_vector_float * tau_V,
                                    gsl_matrix_float * V,
                                    gsl_vector_float * diag, 
                                    gsl_vector_float * superdiag);

int gsl_linalg_bidiag_unpack2_float (gsl_matrix_float * A, 
                                     gsl_vector_float * tau_U, 
                                     gsl_vector_float * tau_V,
                                     gsl_matrix_float * V);

int gsl_linalg_bidiag_unpack_B_float (const gsl_matrix_float * A, 
                                      gsl_vector_float * diag, 
                                      gsl_vector_float * superdiag);

/* Balancing */

int
gsl_linalg_balance_columns_float (gsl_matrix_float * A, gsl_vector_float * D);


__END_DECLS

#endif /* __GSL_LINALG_H__ */
